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ABSTRACT 

Hydrodynamical models of colliding hypersonic flows are presented which explore the 
dependence of the resulting dynamics and the characteristics of the derived X-ray 
emission on numerical conduction and viscosity. For the purpose of our investigation 
we present models of colliding flow with plane-parallel and cylindrical divergence. 
Numerical conduction causes erroneous heating of gas across the contact discontinuity 
which has implications for the rate at which the gas cools. We flnd that the dynamics 
of the shocked gas and the resulting X-ray emission are strongly dependent on the 
contrast in the density and temperature either side of the contact discontinuity, these 
effects being strongest where the postshock gas of one flow behaves quasi-adiabatically 
while the postshock gas of the other flow is strongly radiative. 

Introducing additional numerical viscosity into the simulations has the effect of 
damping the growth of instabilities, which in some cases act to increase the volume of 
shocked gas and can re-heat gas via sub-shocks as it flows downstream. The resulting 
reduction in the surface area between adjacent flows, and therefore of the amount of 
numerical conduction, leads to a commensurate reduction in spurious X-ray emission, 
though the dynamics of the collision are compromised. 

The simulation resolution also affects the degree of numerical conduction. A flner 
resolution better resolves the interfaces of high density and temperature contrast and 
although numerical conduction still exists the volume of affected gas is considerably 
reduced. However, since it is not always practical to increase the resolution, it is 
imperative that the degree of numerical conduction is understood so that inaccurate 
interpretations can be avoided. This work has implications for the dynamics and emis- 
sion from astrophysical phenomena which involve high Mach number shocks. 

Key words: hydrodynamics - methods: numerical - conduction - Shock waves - X- 
rays:general - ISMijets and outflows - stars:outflows 



1 INTRODUCTION 

Colliding hypersonic flows occur in a number of astrophysi- 
cal environments and over a wide range of scales, e.g. mas- 
sive youn g stellar ob jects [Parkin et al. 2009b), astrophys- 
IS hang et al. 2006; Bonito et al. 2007; 

1 I2OQ7I), colliding wi nd binary systems 
1992l:IPittard"2009'). wind-blown bub- 
bles around evolved stars (see Arthur 200 7,, and references 
there - in), SNe (e.g. iTenorio-Tagle et al.l Il99ll : [DwarkadasI 
l20Q7l\ and the cum u lative outflows from young star clusters 



ical iets (lFallelll99ll 



[Sutherland Bickne 



(CWBs, IStevens et a: 



20081: Rodriguez- Gonzalez et al.l 20081; iReves-Iturbide et al 



(ICanto et al.1 I2OQOI: [Rockefeller et al.1 | 2005l: IWunsch et al 

08l;lF 



2OO9I ) and starburst galaxies ([Strickland StevensI I2OOOI : 



E-mail:parkin@astro. ulg.ac.be 



iTenorio-Tagle et al.|[2003l : iTang et al.ll20Q9l ). Flow collisions 
can be subject to turbulent motions, the growth of linear 
and non-linear instabilities in boundary layers, and in some 
cases a globa l instability of the shocked gas (i.e. radiative 
overstability, IChevalier k Imamurall 19821 ). The combination 
of these effects leads to complex scenarios for which numer- 
ical hydrodynamics has proved to be a useful investigatory 
tool. 

However, in the discretization of the governing equa- 
tions of hydrodynamics, additional terms are introduced 
which are purely numerical in origin. Depending on the or- 
der of the scheme, the appearance of these terms acts to 
disperse or dissipate the solution, and therefore terms such 
as "numerical dispersion", "numerical diffusion" or "artifi- 
cial viscosity" are often used to describe them. The undesir- 
able effects of numerical diffusion are minimized as one uses 
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higher order schemes, though all schemes are only first order 
accurate near discontinuites such as shocks (where flow vari- 
ables as well as the perpendicular velocity component, Vp, 
are discontinuous) and contact discontinuities (where there 
is a density and/or temperature jump but Vp is unchanged). 
Contact discontinuities, and interfaces between different flu- 
ids, create special problems for multi-dimensional hydrody- 
namic codes. Unlike shocks, which contain a self-steepening 
mechanism, contact discontinuities spread diffusively during 
a calculation, and c ontinue to broaden as the calculation 
progresses (see e.g. 'Robertson et al."2005'). Some schemes 
employ an algorithm known as a contact discontinuit y steep- 
ener to limit this diffusion (e.g. iFrvxell et ahlbOOOh . How- 
ever, their use remains controversial, since the algorithm is 
based on empirical values with no physical or mathematical 
basis, and requires some care, since under certain circum- 
stances it can produce incorrect results (i.e. "staircasing" , 
Blondin, private communication). 

Purely numerical effects are most prevalent when there 
are large density and temperature contrasts. Unfortunately, 
these frequently occur in practice, as when radiative cooling 
is effective cold dense regions of gas can form. Such regions 
are also inherently unstable, and compressed interface layers 
may be fragmented resulting in cold dense clumps/filaments 
residing next to hot tenous gas. When modelling such phe- 
nomena, the numerical transfer of heat from hot to cold cells 
can change the behaviour of the shocked gas. In particular, 
hot cells on one side of the contact discontinuity can re- 
duce the net cooling rate of denser gas in adjacent cells on 
the other side of the contact discontinuity, and vice-versa. 
A further concern comes when one derives synthetic emis- 
sion from the simulation output. For instance, due to the 
dependence of thermal emission, artificial heating caused 
by numerical conduction can cause dramatic differences in 
the spectral hardness and the magnitude of the integrated 
luminosity. 

The goal of this work is to provide both a qualitative 
and quantitative analysis of the effects of numerical conduc- 
tion and viscosity on the dynamics and observables from 
colliding flows as a function of the density and temperature 
constrast between the post shock gas. For the purposes of our 
investigation we have performed hydrodynamic simulations 
of colliding flows in plane-parallel and cylindrical geometries. 
In both scenarios the influence of efficient radiative cooling 
and powerful instabilities cause cold dense layers/clumps to 
reside next to hot rarefied gas. We show that the calculated 
X-ray emission from the postshock gas is strongly dependent 
on the parameters of the opposing flows. 

The remainder of this paper is structured as follows: in 
§ [2] we give a description of the hydrodynamics code and 
details of the X-ray emission calculations. In §[3] we present 
model descriptions and results, in § [4] a discussion, and we 
close with conclusions in § [S] 



2 NUMERICAL METHOD 
2.1 The hydrodynamics code 

The collision of hypersonic flows is modelled by numerically 
solving the time-dependent equations of Eulerian hydrody- 
namics in a 2D cartesian coordinate system. The relevant 



equations for mass, momentum, and energy conservation 
are: 

l+V.pu = 0, (1) 

^+Vpuu + VP = 0, (2) 

^ + V.[(,^ + P)u] = (4)^(T). (3) 

Here E — e + ^|u|^, is the total gas energy, e is the internal 
energy, u is the gas velocity, p is the mass density, P is 
the pressure, T is the temperature, and mn is the mass 
of hydrogen. We use the ideal gas equation of state, P = 
(7 — l)/pe, where the adiabatic index 7 = 5/3, to close the 
equations. 

The radiative co oling is c alculated f rom the ME KAL 
thermal plasma code (|Kaastra|[l99 2: Me we et al.lll995l ) dis- 
tributed in XSPEC (vl 1.2.0). The temperature of the pre- 
shock flows is assumed to be maintained at ~ 10^ K through 
photoionization heating. Gas in the shocked region which 
rapidly cools is prevented from cooling below this tempera- 
ture. 

The simulations presented in this work were per- 
forme d using the FLA SH version 3.1.1 hydrodynamics 
code (iFrvxeh et al.ll2QQ(t) which uses the piecewise-parabolic 
method IColella Sz Woodward! (|l984h to solve the hydrody- 
namic equations. This code has been des igned to operate 
with a block-structured AMR grid (e.g. iBerger Oliger 
Il989l) using the PARAMESH package (|MacNeice et al. 
2000) under the message-passing interface (MPI) architec- 
ture. In the models presented in this work the AMR uses 
square blocks consisting of 8^ cells. Details relating to the 
adopted resolution and size of the simulation domain for 
the colliding (plane-parallel) laminar flow and cylindrically 
diverging colliding wind binary models are given in §§ 13.11 
and 13. 2[ respectively. Contact discontinuity steepening was 
used in all calculations, with the standard pa rameter (771 = 
20,7/2 = 0.05, e = 0.01, iFrvxeh et aDl2000h . A customized 
unit has been implemented into the FLASH code for op- 
tically thi n radiative cooling using the explicit method de- 
scribed in IStrickland BlondinI (|l995h . An advected scalar 
variable is included in the hydrodynamics calculations to 
distinguish between the flows. 

2.2 X-ray emission 

To calculate the intrinsic X-ray emission from the simulation 
we assume solar abundances and use emissivities for opti- 
cally thin gas in collisional ionization equilibrium obtained 
from look-up tables calculated from the MEKAL plasma 
code containing 200 logarithmically spaced energy bins in 
the range 0.1-10 keV, and 101 logarithmically spaced tem- 
perature bins in the range 10^ - 10^ K. The advected scalar 
is used to separate the X-ray emission contributions from 
each 

^ Further tests in which 771,772, and e were varied revealed that 
within the limit that the dynamics of the simulation remained 
reasonably unaffected, the numerical conduction effects that we 
discuss were not significantly reduced. 

^ Further tests in which a cut is placed on the fluid dye variable 
when calculating the X-ray emission from each fluid component 
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3 RESULTS 

We have performed two sets of simulations to examine the 
effect of numerical conduction on the gas dynamics and the 
derived X-ray characteristics of colliding flows. The first sce- 
nario is that of colliding plane-parallel flow and the second is 
of the wind- wind collision in a massive stellar binary system. 

In many of the models presented in this work the post- 
shock gas has a temperature which places it near a lo- 
cal minimum in the cooling function, A(T), at a value of 
^ -j^Q-23 g-i ^]^jg ^g^^ utilised to determine 

approximate values for the cooling time, tcooi, and the cool- 
ing length, /cool, of postshock gas: 



3kT 
8nA(T) 



10-^°^ s. 



-2vi 



(4) 



(5) 



Zcooi — V • tcool ~ 10 — cm, 

p 

where T is the postshock gas temperature inK, vg is the gas 
velocity in units of 10^ cm s~^, n is the gas number density 
in cm~^, and p the gas density in g cm~^. Note that n, 
and p are the values immediately preshock. 

In models of CWBs the geometry of the colliding flows 
allows one to define an escape time for postshock flow leaving 
the system, tesc- The strength of cooling in these systems can 
then be characterised by a dimensionless cooling parameter 
which is the ratio of the cooling time to the postshock flow 
time (IStevens et al.|[T992. ). 



tcool _ ^8<^12 
^esc M—Y 



(6) 



where vg is the gas velocity in units of 10^ cm s~^, di2 is 
the binary separation in units of 10^^ cm, and M-r is the 
stellar wind mass- loss rate in units oflO^^Moyr"^. Values 
of X ^ 1 are representative of adiabatic gas, whereas x ^ 1 
indicates that the postshock gas is radiative and will cool 
rapidly. In the following we use the subscripts 1 and 2 to 
indicate the cooling parameter for the postshock wind 1 and 
2 material, respectively. 

Using the Rankine-Hugoniot shock jump conditions we 
can estimate the postshock gas temperature, T, and a cor- 
responding energy for (thermal) X-ray emission from the 
postshock gas as, 

kT l.Uvl keV. (7) 



3.1 Colliding plane-parallel flows 

We introduce the effects of numerical conduction with a 
model of colliding plane-parallel flows, relevant to a wide va- 
riety of astrophysical phenomena. Each flow is of constant 
density and velocity, and they are separated by a vertical 
discontinuity through the centre of the grid. The simulation 
domain extends to Xmax = 2/max = 3 x 10^^ cm and is cov- 
ered with a coarse grid consisting of x x y = 15x 15 blocks, 
with each block consisting of 8 x 8 cells. We allow for 3 



did not prove successful in significantly reducing the contamina- 
tion of the results by numerical conduction. The main trend from 
using a cut was that the softer of the two emission components 
becomes harder and vice- versa. 



nested levels of refinement which gives a minimum cell size of 
3.125 X 10^° cm and an effective resolution of 960 x 960 cells. 
The left and right hand boundaries are inflow conditions 
for flow 1 and 2, respectively, while the upper and lower 
boundaries employ zero gradients. The flows are hypersonic, 
so their ram pressure dominates over their thermal pres- 
sure. The simulations are allowed to run for 10^ s which 
is sufficiently long to allow the postshock gas to cool and 
for the initial conditions to flow off the grid. We note 
that previous models of colliding plane-p arallel flows typ- 



ically focus on the s upersonic regi me (e.g. [ Wa lder & F olin: 



1998 2000l:lFolini Walder 20061: iHeitsch et al 2006 . 200i 



bhnil 



iNiklaus et al.ll2009l : iBaneriee et al.ll2009l ) which is more rel- 
evant to star formation. The region of parameter space 
examined here explores the hypersonic regime, a situation 
which is more relevant to systems such as CWBs, SNe, clus- 
ter/galactic winds, and jets. The simulation results highlight 
the turbulent nature of hypersonic flow collisions, which de- 
serve a detailed investigation. However, for the purposes of 
this work we provide a brief description of the dynamics and 
defer a more indepth analysis to a later paper. 

We have run a series of models where the densities and 
velocities of the flows have been varied while maintaining 
equal ram-pressure for the flows (Table [1]) and keeping the 
resolution constant. Through models SLAB- A to SLAB-C 
we decrease (increase) the preshock gas velocity (density) 
of flow 2. Following this, in models SLAB-C to SLAB-E the 
preshock parameters of flow 2 are kept fixed while the veloc- 
ity (density) of flow 1 is decreased (increased). As such, in 
model SLAB-A the postshock gas of both flows is adiabatic, 
whereas in model SLAB-E the postshock gas is strongly ra- 
diative. The preshock temperature of the flows is 10^ K, 
giving Mach numbers of 67, 134, and 300 for preshock ve- 
locities of 1000, 2000, and 3000 km s"\ respectively. The 
postshock gas is allowed to cool back to 10^ K. In the fol- 
lowing, we discuss the effect of these variations on the gas 
dynamics, and then present the results of X-ray calculations 
which provide a quantitative analysis of the effect of numer- 
ical conduction on observable quantities. 

Model SLAB-A consists of two hypersonic flows with 
identical preshock parameters (vi = V2 = 3 x 10^ cm s~^. 



Pi 



P2 = 1.1 X 10"^"^ g cm-^) 



As the flows collide shocks 



are generated which heat gas to T > 10 K. This gas 
slowly cools so that at later times a thin dense shell of 
cold (T ^ 10^ K) gas is formed near the contact discon- 
tinuity separating the flows. This is Ray leigh- Taylor (RT) 
unstable, and flngers of material from the dense shell be- 
gin to protrude into the hot tenous gas (Fig. [1] top pan- 
els). At early times in the simulation the shock front of 
both flows oscillates dramatically (th i s is the radiative over- 
stabili ty, e.g. [c hevalier & Ima mural 19821: Imamura et al 



"1984'; IWalder Folinil 119961 : iPittard et al.l I2OO5I : iMignone 
i2005|). By the time of the snapshot shown in Fig. [1] these 
oscillations have been damped by waves which reflect back 
from the central cold dense layer out of phase with the orig- 
inal shock wave, disrupting the coherence between waves in 
the oscillating shock, and causing it to deteriorate in to an 
approximate steady state ([Strickland Blondin|[T995h . This 
effect is e nhanced by the distorti on of the central cold dense 
gas layer (IWalder Folinil 1 19961 ). 

In models SLAB-B and SLAB-C the velocity (density) 
of flow 2 is decreased (increased) while the ram pressure of 
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Figure 1. Simulation snapshots showing density (left column), temperature (middle column), and flow 2 temperature (right column). 
Models shown from top to bottom: SLAB-A, SLAB-B, SLAB-C, SLAB-D, and SLAB-E. Model parameters are listed in Table [T] Large 
tick marks correspond to a distance of 1 x IQ-*^^ cm. 
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Table 1. Parameters pertaining to the colliding laminar flow simulations. The cooling time, tcoob cind the cooling length, ^coob cire 
calculated from Eqs. |4]and[5l respectively. 



Model pi VI ^cooU ^cooU P2 V2 tcool2 ^cool2 

(g cm~3) (10^ cm s~^) (s) (cm) (g cm~3) (10^ cm s~^) (s) (cm) 



SLAB-A 


1.1 X 10- 


14 


3 


8.18 X 10^ 


2.45 X 10^3 


1.1 X 10- 


14 


3 


8.18 X 10^ 


2.45 X 10^3 


SLAB-B 


1.1 X 10- 


14 


3 


8.18 X 10^ 


2.45 X 10^3 


2.5 X 10- 


14 


2 


1.60 X 10^ 


3.20 X 10^2 


SLAB-C 


1.1 X 10- 


14 


3 


8.18 X 10^ 


2.45 X 10^3 


1.0 X 10- 


13 


1 


1.00 X 10^ 


1.00 X 10^1 


SLAB-D 


2.5 X 10- 


14 


2 


1.60 X 10^ 


3.20 X 10^2 


1.0 X 10- 


13 


1 


1.00 X 10^ 


1.00 X 10^1 


SLAB-E 


1.0 X 10- 


13 


1 


1.00 X 10^ 


1.00 X 10^1 


1.0 X 10- 


13 


1 


1.00 X 10^ 


1.00 X 10^1 




1e+22 I ' • ' ' ' • ' — • — I 

1 10 

Energy (keV) 

Figure 2. X-ray spectra calculated from flow 1 in models SLAB- 
A, SLAB-B, and SLAB-C. Model parameters are noted in Table[T] 
Only flow 2 parameters are varied between the models. 



the preshock flow remains the same. As in model SLAB-A 
the postshock gas of both flows demonstrates over-stability. 
As the preshock velocity (density) of flow 2 decreases (in- 
creases) the frequency of oscillations and the width of the 
layer of postshock gas increase and decrease respectively, as 
tcooi and /cool both decrease. The thin dense shell of gas 
which separates the shocks becomes noticeably more struc- 
tured than in model SLAB-A. At the end of the simulations 
there is a complicated spatial distribution of gas in which 
cold and dense material, lower density shock heated gas, and 
hot, rarefied gas produced by the vigorous action of instabil- 
ities all reside adjacent to one another (Fig. [T]). The hottest 
postshock gas of flow 2 occurs just behind the shock, as is 
also the case in model SLAB-A. 

In model SLAB-C, the postshock gas of flow 2 now cools 
rapidly to form a cold dense shell, and unlike models SLAB- 
A and SLAB-B the width of the region of postshock gas is 
not resolved. Of particular note is the fact that the maxi- 
mum temperature attained by flow 2 gas in model SLAB-C 
is 6.7 X 10^ K, which is substantially higher than the 
expected value of T ^ 1.4 x 10^ K (from the preshock veloc- 
ity). Examining the location of flow 2 gas with T > 10^ K 
(Fig. [1]) we see that the hottest postshock gas now occurs at 
the contact discontinuity between the postshock gas of flows 
1 and 2 in contrast to what is seen in models SLAB-A and 
SLAB-B. This is clearly the effect of numerical conduction 
rearing its ugly head. 

In models SLAB-A to SLAB-C we have examined the 




100000 I ' — ' — ' — 

100000 1e+06 1e+07 1e+08 

Temperature (K) 



Figure 3. Distribution of flow 1 mass as a function of postshock 
gas temperature in models SLAB-A, SLAB-B, and SLAB-C. The 
differences between the distributions at T < 4 x 10^ K are due 
to differing degrees of numerical conduction, whereas at T ^ 4 x 
10^ K they are related to the orientation of the shock with respect 
to the upstream flow. 



effect of progressively increasing the density and tempera- 
ture constrast between the postshock gas whilst keeping the 
parameters of flow 1 fixed so as to have hot, and largely 
adiabatic, postshock gas on at least one side of the contact 
discontinuity. In models SLAB-C to SLAB-E we now look at 
the opposite situation in which we keep highly radiative gas 
on at least one side of the contact discontinuity while reduc- 
ing the preshock velocity of the flow on the other side of the 
contact discontinuity until its postshock gas also becomes 
highly radiative. Through an intermediary model (SLAB- 
D) we arrive at model SLAB-E, in which two equal flows 
with vi = V2 = I X 10^ cm s~^, pi = p2 = 1 x 10~^^ g cm~^ 
produce highly radiative postshock gas. Through this se- 
quence of models the decrease (increase) in the preshock 
velocity (density) of flow 1 causes a corresponding decrease 
in the cooling time and cooling length of its postshock gas. 
Examining progressive trends we note that through models 
SLAB-C to SLAB-E the reduced cooling length in the post- 
shock gas causes the width of the region of postshock gas 
of flow 1 to become narrower and pertubations to the cold 
dense layer increase in ferocity (Fig. [T]). In models SLAB- 
C and SLAB-D the distortion of the cold dens e layer is 
due to a combination of thin-shell (IVishniadllQSah . RT, and 
Kelvin-Helmholtz (KH) instabilities. In model SLAB-E the 
shocks which bound the cold dense layer are isothermal with 
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Figure 4. X-ray spectra calculated from flow 2 in models SLAB- 
C, SLAB-D, and SLAB-E. Model parameters are noted in Table[T] 



large amplitude oscillations resulting from the non-linear 
thin shell instability (NTSI, I Vishniad [Tooi ) . It is interest- 
ing to note that in model SLAB-E almost all of the post- 
shock gas cools rapidly and collapses into the dense layer 
at T < 10^ K. However, there is a small amount of gas at 
T :^ 1.6 X 10'^ K, higher than the 1.36 x lO'^ K expected from 
the flow velocity. We attribute this diflFerence to instability 
driven oscillations moving upstream into the incoming flow, 
and which thus cause higher preshock velocities in the frame 
of the shock. 

To highlight the impact of numerical conduction on the 
inferred observables we examine the 1-10 keV X-ray spec- 
tra calculated from the models (see § 12.21 for details) . In 
Fig. [2] the X-ray spectra calculated from flow 1 in models 
SLAB-A, SLAB-B, and SLAB-C are shown. Comparing the 
spectra from models SLAB-A and SLAB-B we see that as 
V2 is decreased the normalization of the spectrum decreases. 
However, a further reduction in V2 between models SLAB- 
B and SLAB-C actually produces a higher normalization 
in the spectrum of model SLAB-C. To explain these diflFer- 
ences one must examine the temperature plots in Fig. [T] and 
the mass- weighted temperature distributions in Fig. [S] The 
lower normalization in the model SLAB-B spectrum com- 
pared to that of models SLAB-A and SLAB-C is due to a 
lower amount of gas at ^ 10^ K in the former, which is ev- 
idently due to the orientation of the shock with respect to 
the upstream flow - in models SLAB-A and SLAB-C more of 
the shock front is normal to the preshock flow, with more se- 
vere distortions in model SLAB-B. The differences between 
the spectra for models SLAB-A and SLAB-C are the direct 
result of increased numerical conduction in model SLAB-C 
- the higher surface area for interactions in model SLAB-C 
results in a larger amount of energy being extracted from 
hot flow 1 gas by cold flow 2 gas in the temperature range 
T:^ 10^ - 4 X 10^ K. 

We now turn our attention to the effect that numerical 
conduction has on the spectrum of the colder, denser post- 
shock gas. In models SLAB-C, SLAB-D, and SLAB-E the 
preshock parameters of flow 2 are identical. Therefore, intu- 
itively one would expect that the spectrum calculated from 
this gas should also be identical. Yet, examining Fig. [J] we 
see that this is clearly not the case. Comparing the spectra 
calculated for models SLAB-C and SLAB-E, the numerical 



conduction of heat in the former results in a higher normal- 
ization to the spectrum 2 dex). This represents a cause 
for concern as it is clear from the flow 2 temperature plot 
in Fig. [1] that the volume of gas contaminated by numeri- 
cal conduction in model SLAB-C is relatively small, yet the 
impact on the derived X-ray spectrum is significant. 



3.2 Colliding winds in a binary star system 

Colliding stellar winds in binary systems are particularly 
useful for the study of the effect of instabilities on the global 
fiow dynamics and resulting emission characteristics. This is 
in part due to the inherent geometry of the fiows; close to 
the line-of-centres between the stars the stellar winds collide 
head on, postshock gas at the stagnation point is then accel- 
erated in the tangential direction, and there can be consid- 
erable shear between the fiows either side of the contact dis- 
continuity. In addition, the importance of radiative cooling 
in some systems can bring about density and temperature 
contrasts of many orders of magnitude between adjacent gas. 

The model consists of two stars with hypersonic, 
isotropic stellar winds which collide at their terminal speeds. 
We restrict our inv estigation to 2D and neglect any orbital 
motion effects (e.g. iLemaster et al.l l2007l : IParkin PittardI 
l2QQ8l : |Pittardll2QQ9l ). Normally, 2D simulations of colliding 
stellar winds assume axis- symmetry, and an rz grid is used. 
However, such simulations ca n be susceptible to the "car- 
buncle" instability (see e.g. Pittard et al.l 1 19981 ) which is 
purely numerical in origin (|Quirklll994 V To avoid this, we 
have instead adopted slab-symmetry, in which the diver- 
gence of the fiow goes as r~^, and placed the 'stars' in the 
centre of the grid. The grid extends to x = ±1.5 x 10^^ cm, 
y = 0-2.25 X 10^^ cm with star 1 centered at (0, 5x 10^^ cm) 
and star 2 centered at (0, 1 x 10^^ cm), giving a binary 
separation, dsep = 5 x 10^^ cm. For the models presented 
in §§ 13.2.11 and 13.2.21 the grid is initialized with x x y — 
32 X 24 blocks. We allow for 2 nested levels of refinement, 
where adjacent levels differ in resolution by a factor of two. 
This provides an effective resolution on the finest grid of 
1024 X 768 cells. The grid resolution is varied in § [3T3l 
where an investigation into resolution effects is performed. 

To initiate the stellar winds appropriate hydrodynamic 
variables (i.e. />, P, v) are mapped into cells residing within 
a radial distance of the respective star of 5 x 10^^ cm. Model 
parameters are noted in Table [2] The unshocked winds are 
initialized with T — 10^ K. Noting that in 2D cartesian ge- 
ometry the divergence of the fiow goes as we calculate 
the stellar wind density profiles using the following equa- 
tions. 



P2 



Mia 
27TriVi ' 

M2a 



27Tr2V2y/fi 

where 

a — — —-^ — . 



(8) 



(9) 



(10) 



r and v are the radial distance from, and stellar wind veloc- 
ity adopted for, the respective star, and 77 (= M2V2 / Mivi) is 
the wind momentum ratio which is kept constant at a value 
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Figure 5. Simulation snapshots showing density (left column), temperature (middle column), and wind 1 temperature (right column). 
Models shown from top to bottom: CWB-A, CWB-B, CWB-C, CWB-Avisc, and CWB-A+. Wind 1 is flowing from the bottom left 
corner, while wind 2 occupies the top half of each plot. These calculations were performed in slab-symmetry (see text for details). Model 
parameters are listed in Tabl£^_Large tick marks correspond to a distance of 5 x 10-'^'^ cm. For illustrative purposes, only a subsection 
computa^i^na^grKpis' Bicrmi, with the line of centres between the stars actually running through the middle of the grid. 
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of 0.2 in all models. In the above equations we are effectively 
fixing the stagnation point wind density in our 2D cartesian 
geometry simulations to the equivalent attained from a 2D 
axis- symmetric or 3D simulation. This approach has the ad- 
vantage that the preshock gas densities only differ slightly 
from those of an equivalent simulation with an diver- 
gence (i.e. 3D, or 2D axis-symmetric), and as such radiative 
cooling is similar. 

For our fiducial m odel we choose param eters similar 
to those determined by I Parkin et al.l (|2Q09al ) for the mas- 
sive star binary system 77 Car. In the context of CWBs, 
T] Car presents us with some unique problems when it comes 
to modelling its wind-wind collision which we do not find 
when modelling other long-period binaries (e.g. WR140, 
L Ori) . These problems arise from the fact that the slow and 
dense primary wind, once shocked, cools rapidly into a dense 
sheet, whil e the shocked secondary wind behaves largely ad i- 
abatically (IPittard CorcoranI [2QQi : IParkin et al.1 l2Q09al \ 
which produces large jumps in the flow variables across the 
contact discontinuity, with cold dense gas adjacent to hot 
rarefied gas. Inevitably there is some numerical conduction 
of heat across the interfacjfl- 

In the following sections we present results for simula- 
tions performed with low and high numerical viscosity, and 
also the results of resolution tests. In each instance we first 
discuss the gas dynamics and then present the results from 
X-ray calculations. 



3.2.1 Low numerical viscosity 

In model CWB-A wind 1 is highly radiative (xi = 0.01, 
see Eq. (6]) and wind 2 is largely adiabatic (x2 — 400) . The 
postshock wind 1 gas cools rapidly and collapses to form a 
thin dense shell (Fig. [5]). The rapid onset of RT, KH, and 
thin-shell instabilities fragments this shell, causing clumps 
of wind 1 material to become interspersed with the hot post- 
shock wind 2 material. The fragmentation of the interface 
changes the surface area between hot and cold gas, and in- 
creases the total amount of numerical conduction in the cal- 
culation. Postshock wind 2 material reaches peak tempera- 
tures of ^ 10^ K, consistent with the value expected from 
Eq.[7l However, wind 1 gas reaches temperatures of ^ 10^ K, 
much higher than the expected T ^ 3x 10^ K (see the wind 1 
temperature plot in Fig. [5])- As can be seen, the postshock 
region is turbulent, which complicates deciphering heating 
through flow collisions and spurious heating brought about 
by numerical conduction. In ^ I3.2.2l we elucidate the relative 
contributions of these two mechanisms by performing tests 
with additional numerical viscosity to suppress the growth 
of instabilities. 

The intrinsic X-ray spectrum from model CWB-A is 
dominated by emission from the shocked wind 1 material 
at ^ < 2 keV and by wind 2 material above this energy 
(Fig.[6l). The difference in the spectral slopes does not read- 
ily relate to the contrast in preshock velocities. The preshock 



^ To circumvent this problem sophisticated models of CWBs have 
been developed which, at the cost of providing limited information 
abo ut the flow dynamics, avoid details of mixing in the postshock 
gas (lAntokhin et al.ll2004l : [Parkin Pittardll2008r) . 






Figure 6. X-ray spectra calculated from the OWE simulations 
noted in Table [2] The total (solid red line), wind 1 (dashed green 
line), and wind 2 (dotted blue line) emission contributions are 
shown. From top to bottom: CWB-A, CWB-B, CWB-C, CWB- 
Avisc, and CWB-A+. Corresponding integrated luminosities are 
listed in Table E] 
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Table 2. Parameters pertaining to the CWB simulations of ^ 13.21 The wind momentum ratio, 77, is kept constant in all simulations. 



Model 


Comment 


Ml 


VI 


XI 


M2 


V2 


X2 






(lO-^Moyr-i) 


(108 cm s-i) 




(lO-^Moyr-i) 


(108 cm s-i) 




CWB-A 




3000 


0.5 


0.01 


100 


3 


405 


CWB-B 


X2 i 


3000 


0.5 


0.01 


300 


1 


1.67 


CWB-C 


XI t 


500 


3 


81 


100 


3 


405 



velocity of wind 1 (vi = 500 kms~^) corresponds to an en- 
ergy of ^ 0.3 keV (using Eq. (T)). Therefore, we would not 
expect emission extending to energies of ^ 10 keV at the 
magnitude observed in Fig. (6] The total luminosity is domi- 
nated by postshock wind 1 gas (see Table [3|) which produces 
considerable soft X-ray emission. 

Lowering the velocity of wind 2 to = 1 x 10^ cm s~^ 
(model CWB-B), reduces the cooling parameter to a value 
which corresponds to gas which is highly radiative (x2 = 
1.67). Examining Fig. [5] we see that such a change causes a 
dramatic difference to the flow dynamics. The structure of 
postshock gas is now dominat ed by lar^e amplitude oscilla- 
tions caused by the NTSI (e.g. IStevens et al ] |l992l : IVishn"iI3 
I1994I ). Heat conduction from wind 2 into wind 1 is less- 
ened in model CWB-B compared to model CWB-A due to 
a lower temperature gradient at the interface between the 
winds and this has a catastrophic effect on the resulting 
emission (Fig. [6|). The spectrum and total luminosity are 
now dominated by emission from wind 2, and there is a sub- 
stantial reduction in the wind 1 emission at all energies. This 
is unsurprising when one compares the wind 1 temperature 
plots for models CWB-A and CWB-B; for model CWB-B 
there is very little wind 1 gas at T > 10^ K. 

We have explored above the effect that the thermal 
properties (e.g. temperature) of the wind 2 postshock gas 
has on the resulting emission from wind 1, but what effect 
does wind 1 have on wind 2 gas? In model CWB-C the pa- 
rameters of wind 1 have been modified so as to produce 
an adiabatic postshock flow (xi = 81). Compared to mod- 
els CWB-A and CWB-B the WCR now appears relatively 
smooth with only small perturbations of the contact discon- 
tinuity due to velocity shear between the postshock winds 
driving a KH instability (Fig. [5|). Examining the spectra, 
we see that the slope of both the wind 1 and 2 spectra are 
now identical, as expected for identical preshock velocities. 
A small increase in the normalization of the spectrum and in- 
tegrated luminosity calculated from wind 2 is seen between 
models CWB-C and CWB-A, consistent with heat being 
conducted out of wind 2 by the neighbouring cold gas of 
wind 1 in model CWB-A. 

In summary, models CWB-A, CWB-B, and CWB-C 
demonstrate that differences in the density and postshock 
cooling rate across the contact discontinuity are important 
for the dynamics and the calculated X-ray emission. To ex- 
amine how the contamination of the X-ray emission varies 
with the wind 2 parameters we have performed further tests 
which are shown in Fig. [Tj As the preshock velocity of 
wind 2, V2, and therefore X2, is increased the integrated lu- 
minosity of wind 1 and 2 in the hard band (7-10 keV) also 
increases. There is a clear jump in the hard band luminos- 
ity which occurs between ^2 ^ 1 and 1.5 x 10^ cm s~^, and 



□ 

□ 

+ 



□ 

o 



>K >K >K 
□ □ □ 

6 i * 
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Wind 2 
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Wind 2 - viscous Q 
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Wind 2 velocity (10^ cm s'^) 



3.5 



Figure 7. Hard band (7-10 keV) luminosity calculated from 
wind 1 and 2 as a function of V2 for models CWB-A and CWB- 
Avisc. For these model calculations the wind 1 parameters were 
set to those of model CWB-A (-ui = 5 x lO'^ cm s"-*^) and wind 2 
parameters were varied accordingly while maintaining a constant 
wind momentum ratio. In the ideal situation where there is no 
numerical heat conduction between wind 1 and wind 2, the hard 
band luminosity of wind 1 would be negligible. 



Table 3. Integrated 1-10 keV X-ray luminosities calculated for 
the simulated models. Ltot? Li, and L2 are the total, wind 1, 
and wind 2 luminosities (in units of erg s~^ cm~^), respectively. 
Models CWB-Avisc and CWB-A+ are variants of model CWB-A 
with higher numerical viscosity (§ I3.2.2|) and higher simulation 
resolution I3.2.3|) . respectively. 



Model 


Ltot 




Li 






L2 






CWB-A 


9.50 


X 


1020 


6.42 


X 


1020 


3.08 


X 


1020 


CWB-B 


2.42 


X 


1020 


4.44 


X 


1019 


1.98 


X 


1020 


CWB-C 


9.52 


X 


1020 


5.86 


X 


1020 


3.66 


X 


1020 


CWB-Avisc 


4.14 


X 


1020 


2.48 


X 


1020 


1.66 


X 


1020 


CWB-A+ 


7.84 


X 


1020 


5.20 


X 


1020 


2.64 


X 


1020 



corresponds to the transition of wind 2 from a radiative to 
an essentially adiabatic postshock flow. When both winds 
are similar in character (i.e. radiative in this case) there ap- 
pears to be little heat conduction between them, while large 
differences in the wind speeds and postshock temperatures 
result in significant numerical heat conduction, as shown by 
the rise in the hard band luminosity from wind 1 as V2 in- 
creases. 



© 2009 RAS, MNRAS OOO.fTIfM] 



10 E. R. Parkin & J. M. Pittard 



3.2.2 High artificial viscosity 

Hydrodynamical codes typically allow the user to modify 
the magnitude of additional diffusion- like terms which act 
as artificial viscosity in the governing flow equations. Intu- 
itively, one would expect the addition of artiflcial viscosity 
to increase the level of heat conduction across a unit area 
of the interface between cold and hot gas. However, viscos- 
ity also suppresses instabilities, and thus reduces the total 
area for interactions. Whether the total amount of numerical 
heat conduction increases or decreases with increasing arti- 
ficial viscosity will depend on the relative strength of these 
effects. 

Low numerical viscosity in model CWB-A results in the 
growth of instabilities and thus a turbulent wind-wind col- 
lision region (WCR). The inclusion of additional numerical 
viscosity in model CWB-Avisc surpresses the growth of such 
instabilities, and in comparison to model CWB-A the WCR 
is much smoother (Fig. [5]). The volume occupied by shocked 
wind 2 material is significantly smaller in model CWB-Avisc 
than in model CWB-A. There appears to be two reasons for 
this. The reduction in volume far downstream of the apex 
of the WCR is caused by the lack of secondary shocks in 
this region of the fiow, which in model CWB-A result from 
the intrusion of dense clumps of gas from wind 1 into the 
shocked wind 2 gas. However, there is also a significant re- 
duction in the width of postshock wind 2 gas at the apex 
of the WCR which is evidence for the enhanced numerical 
transport of heat out of the gas by the artificial viscosity. 

A comparison of the X-ray calculations for models 
CWB-A and CWB-Avisc reveals that the latter is of lower 
luminosity and spectral hardness (see Tableland Fig. [6]). 
This can be understood by the fact that in model CWB- 
A the higher level of interaction between postshock wind 
1 and 2 gas (i.e. slowly moving dense clumps being struck 
by higher speed, hot, rarefied fiow) involves collisions which 
heat wind 1 gas to soft X-ray emitting temperatures and 
also act to heat wind 2 gas, whilst in model CWB-Avisc 
heat is conducted out of the hottest gas near the apex by 
the increased artificial viscosity. This point is highlighted 
by a comparison of the distribution of mass as a function of 
temperature for models CWB-A and CWB-Avisc (Fig. [8|). 
For model CWB-Avisc there is almost three orders of mag- 
nitude less mass at T :^ 10^ K compared to model CWB- 
A. In addition, the distribution of wind 2 material in the 
range T = 10^ - 10^ K drops-off more rapidly for model 
CWB-Avisc than for model CWB-A. These results confirm 
that i) the increased surface area for interactions in model 
CWB-A compared to model CWB-Avisc results in collisions 
which heat additional (mainly downstream) wind 2 mate- 
rial to T 10^ K through further shocks while also heating 
wind 1 material to T < 10^ K through friction and enhanced 
heat conduction (though the mass-weighted temperature is 
actually reduced - see below) , and ii) enhanced heat conduc- 
tion near the apex of the WCR reduces the temperature of 
the wind 2 gas in model CWB-Avisc. 

In model CWB-A the instabilities help to re-heat both 
wind 1 and 2 gas as it fiows downstream resulting in re- 
curring spikes in the mass-weighted postshock gas temper- 
ature as one moves to larger distances from the shock apex 
(Fig. [9|). However, despite this additional heating, as the 
gas fiows off the grid the mass- weighted wind 2 tempera- 



ture for model CWB-A is in fact lower than that of model 
CWB-Avisc. This would suggest that although the volume 
of postshock gas in model CWB-Avisc is smaller in size than 
that of model CWB-A, the enhanced mixing at the CD in 
model CWB-A produces a lower mass- weighted temperature 
than the smooth fiow of model CWB-Avisc. This does not, 
however, lead to a harder fiow 2 spectrum from model CWB- 
Avisc, because the hottest wind 2 gas in the models occurs 
closer to the apex of the WCR (see Fig. [6] and Table [S]). At 
the apex of the WCR the temperature of postshock wind 2 
gas is lower in model CWB-Avisc than in model CWB-A, 
consistent with the expected extra numerical conduction in- 
troduced by additional artificial viscosity, which ultimately 
modifies the gas temperature and pressure. 

The wind 1 distributions in Fig. [9] provide further ev- 
idence for enhanced numerical conduction as the wind 1 
temperature in model CWB-Avisc is clearly higher than in 
model CWB-A. However, it should be noted that in the cal- 
culations presented in Fig. [9] a cut of T = 10^ K was made, 
and thus the distributions shown are for mass above this 
limit - in model CWB-Avisc all postshock wind 1 gas at 
T > 10^ K has an average temperature of 4 x 10^ K 
whereas in model CWB-A this value is :^ 7 x 10^ K. Exam- 
ining the location of gas at these mean temperatures we see 
that in model CWB-Avisc it resides in a thin layer at the 
contact discontinuity whereas in model CWB-A the turbu- 
lent WCR causes this gas to outline wind-clump interactions 
as well as highlight the heavily distorted contact discontinu- 
ity (see the wind 1 temperature plots in Fig. [5|. Although 
there is more wind 1 gas at T > 10^ K in model CWB-A as a 
result of vigorous mixing (see Fig. [8]), in model CWB-Avisc 
heating at the CD via numerical conduction actually causes 
the mass-weighted temperatures to be higher, resulting in a 
slightly harder wind 1 spectrum (see Fig. (6]). 

Focusing now on the inferred hard band (7-10 keV) lu- 
minosity, we have performed a series of simulations with 
high artificial viscosity and with the wind 1 parameters 
fixed, but with different parameters for wind 2, the results 
of which are plotted in Fig. [3 At ^ 2.5 x 10^ cm s~^ 
(X2 ^ 163) the wind 1 and wind 2 data points from the 
low and high viscosity calculations correlate well. However, 
at ^ 2 X 10^ cm s~^ (x2 ^ 53) the wind 1 points 
with and without additional viscosity diverge considerably. 
This result implies that the effects of numerical conduc- 
tion are somewhat (v2 ^ 2 x 10^ cm s~^) to significantly 
{v2 ^ 10^ cm s~^) reduced by the application of additional 
artificial viscosity when instabilities that otherwise would 
occur are strongly suppressed, thereby reducing the surface 
area between fiows. However, the cost of this approach is 
an unrealistic description of the dynamics occuring in the 
WCR, which can effect the derived observational character- 
istics (i.e. the spectral hardness of the wind 2 emission from 
wind 2) as a result of limiting the presence of smaller scale 
downstream shocks. 

3.2.3 Simulation resolution 

To assess the effect of the simulation resolution on the dy- 
namics of the wind-wind collision region and the result- 
ing X-ray emission, calculations have been performed with 
model CWB-A parameters and with cell sizes in the range 
(1.17- 11.7) X 10^^ cm. 
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Figure 8. Distribution of mass as a function of postshock gas 
temperature in model CWB-A (upper panel) and model CWB- 
Avisc (lower panel). 



At lower resolution the coarser grid essentially acts like 
viscosity; the wavelength of resolvable instabilities is larger 
and the timescale for the growth of the smallest resolvable 
instabilities is longer. This affects the formation of structure 
on smaller spatial scales as one tends towards lower resolu- 
tion (i.e. larger cell sizes). To the contrary, at higher reso- 
lution the formation of structure is enhanced and to illus- 
trate this point we show a snapshot of the highest resolution 
simulation (model CWB-A+; twice the resolution of model 
CWB-A) in Fig. [5] The dense shell of cold, postshock pri- 
mary wind fragments into more numerous clumps of smaller 
scale which as before pass into the supersonic stream of post- 
shock wind 2 gas, forming bow shocks and tails. Because of 
the smaller size of these clumps the timescale for their de- 
struction is shorter. This may limit the ability of clumps to 
traverse co mpletely out of the WCR, due to orbital motion, 
as shown bv lPittardl (| 20091 ). Despite the improved simulation 
resolution, wind 1 gas is still being heated to temperatures 
well in excess of those expected from the preshock wind ve- 
locity. 

Comparing the integrated 1-10 keV luminosities from 
models CWB-A and CWB-A+ reveals lower values for the 
latter simulation (Table [3]). The explanation for this can be 
found in the spectrum for model CWB-A+ (Fig.[6|) where we 
see that the comparative decrease in the integrated luminosi- 
ties is due to a lower normalization in the wind 1 emission 
at all energies, and a reduction in soft (^ ^ 3 keV) emission 



100000 I ^ 1 1 

5e+14 1e+15 1.5e+15 

Distance (cm) 

Figure 9. Mass- weighted postshock gas temperature as a func- 
tion of distance from the apex of the WCR in model CWB-A 
(upper panel) and model CWB-Avisc (lower panel). When cal- 
culating the distributions a temperature cut of T ^ 10^ K was 
used. 



from wind 2. To further examine this point we have per- 
formed a resolution test, the results of which are presented 
in Fig. [TOl The general trend is for the artificial heating of 
wind 1 (by numerical conduction) to decrease as the resolu- 
tion increases (cell size decreases). Therefore, the increased 
simulation resolution has the effect of reducing the net emis- 
sion calculated from wind 1, and hardening the spectrum 
from wind 2. 

In conclusion, the degree of fragmentation and the size 
of the clumps is dependent on the adopted resolution. Higher 
resolution simulations will create smaller clumps and vice- 
versa. In many (if not most) astrophysical situations, hydro- 
dynamical codes cannot simulate the large Reynold's num- 
ber flows that occur in reality. Thus the instabilities are 
not as small as they should be. The increasing popularity 
of su b-grid turbulence m odels may alleviate this problem 
(e.g. IPittard et al.l [ioool v Finally we note that magnetic 
fields/pressure can limit the compression/density contrast 
of gas in the dense shell, and may set a minimum size of 
clumps, while also suppressing ablation and yielding longer- 
living clumps. 
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Figure 10. Hard band (7-10 keV) luminosity from wind 1 as a 
function of the grid resolution. Model CWB-A parameters were 
used for these calculations (see ^ I3.2.3|) . 



4 DISCUSSION 

The numerical conduction of heat between hot and cold cells 
is an undesirable side-effect of the approach used to solve the 
governing equations of fluid flow. It has implications for a 
wide range of phenomena, examples of which include: 

• Colliding flows: the collision of opposing super-sonic 
flows is a common occurance in scenarios ranging from 
the formation of raolecul a r clouds (iFolini WaldeJ l2006l : 
iHeitsch et aD l2QQ6l . l2008l : iNiklaus et al.1 l2009l\ the inter- 
actions of stellar winds from YSOs ([Pelamarter et al.l 
'2000: 'Parki n et al.ll2QQ9bh . and CWBs (IStevens et al.lll992l : 
.Pittard 200^)" where the preshock w i nds may also be in- 
homogeneous (IWalder Folinil l2002l : IPittardI l2007l . l2009l : 
IPittard Parkinll2010l )." Thermal instability and the cool- 
ing of gas are of key importance to the early stages of star 
formation, and numerical heat conduction will have conse- 
quences for simulation results. 

• Expanding bubbles: the expansion of a high pressure 
bubble which sweeps up colder surrounding material 
into a dense shell is releva nt for wind-b lown bubbles 
around massive stars (see lArthud l2007l . and refer- 
ences t here- in) and you n g stel lar clusters (ICanto et al.l 



20001: 



Rockefeller et al.l 



Rodriguez- Gonzalez et aU 
20091). 



ex p andin g thin-shells 

(iDale et al.1 120091) and SNe (e 
Il99ll : iDwarkadasI l2007l : iFerrand et al, 



IWiinsch et al.1 l2008l: 
iReves-Iturbide et al.l 
aro und HI I regions 
Tenor io-Tagle et al.l 



■_ — ■ ■ , and star- 
burst regions (IStrickland StevensI l2000l: iMarcolini et ahl 



I2OIOI). 



20041: ICooper et al.1 I2OO8I: IStrickland HeckmanI l2009l : 
Fuiita et al.ll2009l : iTang et al.ll2009l ). Numerical conduction 
at the interface between the hotter gas in the interior of the 
bubble, SNR, or HII region, and the surrounding cold dense 
shell, leads to the evaporation of material from the shell 
into the bubble, and will change the emission and dynamics 
of these objects. 

• Flow-clump interactions: the interaction of a fast (and 
often hot) diffuse flow and a clumpy medium is a com- 
mon occurance in as tronomy, occuring in such diverse ob- 
jects as PNe (e.g. iMeaburn et al. 1 Il998l : fnlt suura et al 
2009) and starburst s uperwinds (e.g. [Strickla nd et al 



ulations of this interaction now exists (e.g. Klein et al. 


1994 


: Mac Low et al.1 1994: ICregori et al.1 1999. 


200d; 


Mellema et al.1 12OO2I: iFalle et al.1 12OO2I: iFragile et al. 


2003, 


2005 


: iMelioh & de Gouveia Dal Pinol 20041: iMelioh et al. 


2005 


: lOrlan 


do et al.1 I2OO5. I2OO6I. I2OO8I: IPittard et al. 


2005 


. I2OO9I: 


Dvson et al.1 I2OO6I: iTenor io-Tagle et al.1 20061: 


van Loo et a 


[.1 2007. I2OIO 


: Shin et al. 2008: Pittard et al. 


20091: Cooper et al.ll2009l: 


Yirak et al.l 20091). Observational 



signatures such as Ha, OVI, and X-ray emissio n from star- 
burst regio ns (e.g. IWestmoquette et al.1 l2007l) and super- 
winds (e.g. IStricklandet al.ll2000l : ICecil et al.|[2 QQ2l) will be 
sensitive to the level of heat conduction between the hot and 
cold phases, artificial or otherwise. 

The X-ray spectra calculated from the colliding plane- 
parallel flow models presented in § 13.11 show that large de- 
viations can arise from numerical conduction. For example, 
emission from flow 2 in model SLAB-C is 100 times greater 
than that in model SLAB-E. Numerical heat conduction is 
also seen in the CWB models. Of course, in reality some 
physical mixing/diffusion may occur. In this work we have 
not considered thermal electron conduction which, due to 
the large temperature gradients present in the simulations, 
may be important. For instance, the transport of heat across 
the contact discontinuity by th ermal electrons can "evapo- 
rate" cold dense postshock gas (iMvasnikov Zhekovl [19981 : 
IZhekov M vasnikov 1991), affecting the flow dynamics and 
the derived X-ray emissivities. However, depending on the 
strength and orientation of the magnetic field the efficiency 
of th ermal electron condu ction may be significantly hindered 
(e.g. lOrlando et al.|[2008l ). Therefore, to accurately account 
for the effect of thermal electron conduction requires the in- 
clusion of the relevant electron transport physics and mag- 
netic fields in the simulations. With this in mind, it is im- 
portant also to understand the level at which numerical con- 
duction acts at. 

The current investigation has focused on grid-based 
hydrodynamics (GH), of whi ch there are numer o us codes 
availa ble to the communitv ([Frvxell et al.l 2000l: Norman 
200d: iTevssierl [iooi; lO 'Shea et al.1 I200I : Mignone et"al 



I2OO7I : IStone et al.1 I2OO8I ). This choice is in part justified 



by the finding that GH is considerably better at handling 
stron g shocks, contact discontinuities, and in stabilities than 
SPH (lAgertz et al.ll2007l : ffasker et al.ll2008l ). although re- 
cent developments to the SPH scheme have improved its 



20071: Pried 20081: iRead 


et al. I2OO9I: Kawata et al.1 20091: 


Saitoh k Makinol l2009l: 


Cha et al.1 I2OIOI). In hght of the 



growing literature of code comparisons which aim to elu- 
cidate differences between simulation results from GH 



and SPH codes (e.g. Frenk et al. Il999l: lAgei 


^tz et al. 2007: 


Commercon et al.1 20081: iTasker et al.ll2008: 


Wadslev et al.l 


20081: iKitsionas et al. I2OO8I). it would be verv interesting to 



2OOOI : ICecil et al.1 12OO2I ). A large body of numerical sim- 



repeat the current investigation with a state-of-the-art SPH 
code. 



5 CONCLUSIONS 

We have presented hydrodynamic models of colliding hyper- 
sonic flows with the aim of examining the effects of numeri- 
cal conduction on the simulation dynamics and the derived 
X-ray characteristics. The conduction of heat occurs across 
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flow discontinuities due to diffusive terms introduced in the 
discretization of the governing flow equations. A key con- 
clusion from this work is that the magnitude of the numer- 
ical heat conduction is strongly related to the density (and 
temperature) contrast between adjacent gas. X-ray calcula- 
tions performed on the simulation results show that signifi- 
cant changes to spectra can occur by numerical conduction 
alone. Further tests performed with additional artificial vis- 
cosity reveal a complicated relationship between the flow 
dynamics, the magnitude of numerical conduction, and the 
resulting X-ray emission. For instance, the inherent insta- 
bility of the collision regions of hypersonic flows naturally 
enhances the interface area between the flows, which in turn 
enhances the level of numerical conduction. Introducing suf- 
ficient viscosity to damp the growth of instabilities can re- 
duce these effects, but the additional diffusion introduced 
into the fluid equations may increase the level of numerical 
heat conduction where the interface is relatively stable (e.g. 
near the apex of the wind- wind collision region in a colliding 
winds binary system). Finally, we note that while enhanc- 
ing the resolution of the simulation increases the growth of 
small scale instabilities, and thus the area of the interface 
between the hot and cold phases, the overall effect of nu- 
merical conduction is reduced. 

In the present work we have highlighted a fundamen- 
tal problem encountered when using grid-based hydrody- 
namics to model fluids where high density and temperature 
contrasts are present - conditions which can be found in a 
multitude of astrophysical phenomena. Unfortunately, there 
is no simple fix. The brute-force approach to resolving this 
problem would be to employ higher simulation resolution, 
though this is not always a realistic option. 
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